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from concentrated sources in turbulent channel flow 
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Dispersion of a passive scalar from concentrated sources in fully developed turbulent channel 
flow is studied with the probability density function (PDF) method. The joint PDF of velocity, 
turbulent frequency and scalar concentration is represented by a large number of Lagrangian parti- 
cles. A stochastic near-wall PDF model combines the generalized Langevin model of Haworth and 
Pope^ with Durbin's^ method of elliptic relaxation to provide a mathematically exact treatment of 
convective and viscous transport with a non-local representation of the near-wall Reynolds stress 
anisotropy. The presence of walls is incorporated through the imposition of no-slip and imperme- 
ability conditions on particles without the use of damping or wall-functions. Information on the 
turbulent timescale is supplied by the gamma-distribution model of van Slooten et al.'^. Two dif- 
ferent micromixing models are compared that incorporate the effect of small scale mixing on the 
transported scalar; the widely used interaction by exchange with the mean (lEM) and the interac- 
tion by exchange with the conditional mean (lECM) model. Single-point velocity and concentration 
statistics are compared to direct numerical simulation and experimental data at Rct — 1080 based 
on the friction velocity and the channel half width. The joint model accurately reproduces a wide 
variety of conditional and unconditional statistics in both physical and composition space. 

Keywords: turbulent channel flow; scalar dispersion; probability density function method; particle method; 
Langevin equation; Monte-Carlo method 



I. INTRODUCTION 

In engineering industry and atmospheric transport and 
dispersion modeling there is an increasing use of com- 
putational methods to calculate complex turbulent flow 
fields. Many of these computations depend on the k-e 
turbulence model^'^, while some are based on second- 
moment closures^"^. The aim of these statistical meth- 
ods is to predict the first and second moments of the 
turbulent velocity field, respectively. In large-eddy simu- 
lation (LES) the large-scale three-dimensional unsteady 
motions are represented exactly, while the small-scale 
motions are parameterized. As long as the transport- 
controlling processes of interest (eg. mass, momentum 
and heat transfer in shear flows) are resolved, LES pre- 
dictions can be expected to be insensitive to the de- 
tails of residual-scale modeling. In applications such 
as high-Reynolds-number turbulent combustion or near- 
wall flows, however, where the important rate-controlling 
processes occur below the resolved scales, the residual- 
scale models directly influence the model predictions. 
Since there is no universally 'best' methodology that is 
applicable for every type of practical flow, it is valuable 
to develop improvements for the full range of turbulence 
modeling approaches. 

The development of probability density function 
(PDF) methods is an effort to provide a higher-level sta- 
tistical description of turbulent flows. The mean velocity 
and Reynolds stresses are statistics of (and can be ob- 
tained from) the PDF of velocity. In PDF methods, a 
transport equation is solved directly for the PDF of the 
turbulent velocity field, rather than for its moments as 
in Reynolds stress closures. Therefore, in principle, a 



more complete statistical description can be obtained. 
While for some flows (e.g. homogeneous turbulence) this 
higher-level description may provide little benefit over 
second moment closures, in general the fuller descrip- 
tion is beneficial in allowing more processes to be treated 
exactly and in providing more information that can be 
used in the construction of closure models. Convection, 
for example, can be exactly represented mathematically 
in the PDF framework, eliminating the need for a clo- 
sure assumption^'^. Similarly, defining the joint PDF of 
velocity and species concentrations in a chemically reac- 
tive turbulent flow allows for the treatment of chemical 
reactions without the burden of closure assumptions for 
the highly nonlinear chemical source terms^-'^. This latter 
advantage has been one of the most important incentives 
for the development of PDF methods, since previous at- 
tempts to provide moment closures for this term resulted 
in errors of several orders of magnitude^^ . 

In the case of turbulent flows around complex geome- 
tries the presence of walls requires special treatment, 
since traditional turbulence models are developed for 
high Reynolds numbers and need to be modified in the 
vicinity of walls. This is necessary because the Reynolds 
number approaches zero at the wall, the highest shear 
rate occurs near the wall and the impermeability condi- 
tion on the wall-normal velocity affects the flow up to 
an integral scale from the wall^'^. Possible modifications 
involve damping functions^^"^'' or wall-functions^'^"^^. In 
those turbulent flows where a higher level of statistical 
description is necessary close to walls, adequate repre- 
sentation of the near- wall anisotropy and inhomogeneity 
is crucial. Durbin^ proposed a Reynolds stress closure 
to address these issues. In his model, the all-important 
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process of pressure redistribution is modeled through an 
eUiptic equation by analogy with the Poisson equation, 
which governs the pressure in incompressible flows. This 
represents the non-local effect of the wall on the Reynolds 
stresses through the fluctuating pressure terms. In an ef- 
fort to extend PDF methods to wall-bounded turbulent 
flows, Durbin's elliptic relaxation method has been com- 
bined with the generalized Langevin model^ by Drcebcn 
and Pope^^'^'^. With minor simplifications this approach 
is closely followed throughout the present study to model 
the joint PDF of the turbulent velocity field. 

The dispersion of scalars (e.g. temperature, mass, 
etc.) in turbulent fiows is relevant to a number of 
scientific phenomena including engineering combustion 
and atmospheric dispersion of pollutants. Reviews on 
the subject have been compiled by Shraiman and Sig- 
gia^'' and Warhaft^^. Several experimental studies have 
been carried out in order to better understand the be- 
havior of transported scalars in homogeneous isotropic 
turbulence^^"^'^ . A literature review of dispersion from a 
concentrated source in homogeneous but anisotropic tur- 
bulent shear flows is given by Karnik and Tavoularis^^. 
Inhomogeneous turbulence (e.g. the atmospheric bound- 
ary layer or any practical turbulent flow) adds a signif- 
icant level of complexity to these cases. Extensive mea- 
surements of the mean, variance, intermittency, proba- 
bility density functions and spectra of scalar have been 
made by Fackrell and Robins'^° in a turbulent bound- 
ary layer. One point statistics in turbulent channel flow 
have recently been reported by Lavcrtu and Mydlarski'^^, 
whose experimental data are used as the main reference 
point throughout the current numerical modeling study. 

Direct numerical simulation (DNS) has served as an 
important counterpart to measurements of turbulence at 
moderate Reynolds numbers, shedding light on quantities 
that are difficult to measure (e.g. Lagrangian statistics) 
and at locations where nearly impossible to measure (e.g. 
close to walls). Turbulent velocity statistics extracted 
from DNS of channel flow have been reported by Moser 
et al.^^ and Abe et al.^'^, while Vrieling and Nieuwstadt'^^ 
performed a DNS study of dispersion of plumes from sin- 
gle and double line sources. 

A widely used model to incorporate the effects of small 
scale mixing on the scalar in the PDF framework is the 
interaction by exchange with the mean (lEM) modeP'"'''^^. 
While this model has the virtue of being simple and effi- 
cient, it fails to comply with several physical constraints 
and desirable properties of an ideal mixing model^^. Al- 
though a variety of other mixing models have been pro- 
posed to satisfy these properties'^^, the lEM model re- 
mains widely used in practice. Recently, increased at- 
tention has been devoted to the interaction by exchange 
with the conditional mean (lECM) model. Sawford'^* 
has done a comparative study of scalar mixing from line 
sources in homogeneous turbulence employing both the 
lEM and lECM models, wherein he demonstrated that 
the largest differences between the two models occur in 
the near-field. He also investigated the two models in a 



double scalar mixing layer^^ with an emphasis on those 
conditional statistics that frequently require closure as- 
sumptions. Based on the lECM model, PDF micromix- 
ing models have been developed for dispersion of passive 
pollutants in the atmosphere by Luhar and Sawford**^ 
and Cassiani et al.^^"^'^. These authors compute scalar 
statistics in homogeneous turbulence and in neutral, con- 
vectivc and canopy boundary layer by assuming a joint 
PDF for the turbulent velocity field. However, no previ- 
ous studies have been conducted on modeling the joint 
PDF of velocity and a passive scalar from a concentrated 
source in inhomogeneous flows. 

We have developed a complete PDF-IECM model for 
a fully developed, turbulent, long-aspect-ratio channel 
flow, where a passive scalar is continuously released from 
concentrated sources. The joint PDF of velocity, charac- 
teristic turbulent frequency and concentration of a pas- 
sive scalar is computed using stochastic equations. The 
flow is explicitly modeled down to the viscous sublayer 
by imposing only the no-slip and impermeability condi- 
tion on particles without the use of damping, or wall- 
functions. The high-level inhomogencity and anisotropy 
of the Reynolds stress tensor at the wall are captured 
by the elliptic relaxation method. A passive scalar is 
released from a concentrated source at the channel cen- 
terline and in the viscous wall-region. The effect of 
small-scale mixing on the scalar is mainly modeled by 
the lECM model. The performance and accuracy of the 
lECM model compared to the simpler, but more widely 
used lEM model are evaluated. Several one-point uncon- 
ditional and conditional statistics are presented in both 
physical and composition spaces. An emphasis is placed 
on common approximations of those conditional statistics 
that require closure assumptions in concentration-only 
PDF methods, i.e. in methods that assume the under- 
lying turbulent velocity field. The results are compared 
to the DNS data of Abe et al.'^^ and the experimental 
data of Lavertu and Mydlarski^^. The experiments were 
performed at two different Reynolds numbers {Rer = 
Urh/v = 520 and 1080 based on the friction velocity 
Mt-, the channel half width h, and the kinematic viscosity 
v) in a high-aspect-ratio turbulent channel flow, measur- 
ing one point statistics of a scalar (temperature) emitted 
continuously at three different wall-normal source loca- 
tions from concentrated line sources. Measurements were 
performed at six different downstream locations between 
4.0 s$ x/h ^ 22.0. 

The remainder of the paper is organized as follows. In 
Section II, the turbulence and micromixing models are 
described. A brief account of the underlying numerical 
methods, various implementation details together with 
wall-boundary conditions are presented in Section III. In 
Section IV, one-point velocity statistics are compared to 
direct numerical simulation data at Rtr = 1080, and a 
comparative assessment of the two micromixing models 
with analytical and experimental data is also given. De- 
tailed statistics of scalar concentration calculated with 
the lECM micromixing model are presented. Finally, 
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conclusions are summarized in Section V. 



II. GOVERNING EQUATIONS 

The governing equation for the motion of a viscous, 
incompressible fluid is the Navier-Stokes equation 



dt dxj p dxi * ' 



(1) 



where Ui, P, p and u are the Eulerian velocity, pres- 
sure, constant density and kinematic viscosity, respec- 
tively. From Eq. (1) an exact transport equation can be 
derived for the one-point, one-time Eulerian joint PDF 
of velocity f{V;x,t)^°''^-^, 



df 
dt 



dl 

dxi 



Id(P) 



p dxi dVi 
dU, dUj 




U{x,t) = V 



(2) 

where V is the sample space variable of the stochastic 
velocity field U{x,t) and the pressure P has been decom- 
posed into mean (P) and fiuctuating part p. Since the 
PDF f{V; X, t) contains all one-point statistics, mean ve- 
locity and Reynolds stresses are readily available through 



m 

{u,Uj) 



VJdV, (3) 

(v, - ([/,)) (y, - {u,))fdv, (4) 



where the fiuctuation Ui = Vi — (Ui) and the integrals are 
taken over all the three-dimensional sample space of the 
velocity field. While in principle, after modeling the un- 
closed conditional statistics, Eq. (2) could be solved with 
traditional numerical techniques, such as the finite dif- 
ference or finite element method, the high dimensionality 
of the equation prevents an efficient solution with these 
methods. In general, Monte-Carlo techniques are com- 
putationally more efficient for problems with such high 
dimensions. In PDF methods, therefore, a Lagrangian 
description has been preferred. The Navier-Stokes equa- 
tion (1) is widely used to model incompressible fiows in 
the Eulerian framework. An equivalent model in the La- 
grangian framework can be written as a system of gov- 
erning equations for Lagrangian particle locations Xi and 
velocities Ui^'^ 



dX, ^ U,dt + {2iy)^^'^ dW, (5) 
+ 2-^dt+(2.)^/^|^dW(^) 



1 dP 

dUi{t) = — ^dt 
p dx. 



dxjdxj 



dxi 



where the isotropic Wiener process^^ dWi, which is a 
known stochastic process with zero mean and variance 
dt, is identical in both equations (the same exact series 
of Gaussian random numbers) and it is understood that 
the Eulerian fields on the right hand side are evaluated 
at the particle locations Xi. Equation (5) governs the po- 
sition of a stochastic fluid particle which undergoes both 
convectivc and molecular motion. In other words, besides 
convection the particle diffuses in physical space with co- 
efficient h', thus it carries momentum as molecules do 
with identical statistics, as in Brownian motion^^. Each 
of the above three descriptions, i.e. Eqs. (1), (2) and (5- 
6), represents the viscous stress exactly. A remarkable 
feature of the PDF formulation is that the effects of con- 
vection and viscous diffusion, fundamental processes in 
near- wall turbulent flows, have exact mathematical rep- 
resentations, therefore need no closure assumptions. If 
Reynolds decomposition is applied to Eq. (6) 



dt + zv— — - — dt 



p dx 

-^dt + 2v 

p dXi 



dxjdxj 



(2^) 



1/2 



dxi 



dWi 



dxj dxj 



dui 



dxj 



(7) 



the last three terms are unclosed. Closure hypotheses for 
these terms can be made either by suggesting approxima- 
tions for the PDF fluxes (represented by the conditional 
expectations) in Eq. (2) or by proposing stochastic pro- 
cesses that simulate the physical phenomena represented 
by the second row of Eq. (7). For an overview of the wide 
variety of modeling approaches see the review compiled 
by Dopazo'^''. To model the increments in particle veloc- 
ity the generalized Langevin model (GLM) of Haworth 
and Pope^ is adopted here 



(2^)' 



dUi[t) ^ — dt -I- 2i/- — - — dt 

p dxi dXjdXj 

+ {Uj - {U,))dt + [Cae f^ dWi 



/2 dm 



dWi 



(8) 



where Gij is a second-order tensor function of velocity 
statistics, Cq is a positive constant, e denotes the rate 
of dissipation of turbulent kinetic energy and dW- is an- 
other Wiener process. By comparing Eqs. (7) and (8), it 
is apparent that the terms in Gij and Co jointly model 
the last three terms in Eq. (7) representing pressure redis- 
tribution and anisotropic dissipation of turbulent kinetic 
energy A particular specification of Gij corresponds to 
a particular Lagrangian stochastic model for the instan- 
taneous particle velocity increment. Following Dreeben 
and Pope^'^, we specify Gij and Co employing Durbin's 
elliptic relaxation method using the constraint 



{1 + ^Co) e + Gij (u^uj) =0, 



(9) 



which ensures that the kinetic energy evolves correctly in 
homogeneous turbulence^'^. Introducing the tensor pij to 
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characterize the non-local effects G;,- and Cn are defined 



as 



Mil and Co = ^%^, (10) 
: 6ke 



where fc = i {uiUi) represents the turbulent kinetic en- 
ergy. The non-local term pij is specified with the follow- 
ing elliptic relaxation equation 



k (w) % + kH,jki- 



djUk 
dxi 



where 



and 



UkUk 



(11) 



(12) 



(13) 



(14) 



is the Reynolds stress anisotropy, (w) denotes the mean 
characteristic turbulent frequency and Ci, 6*2,75, C„ are 
model constants. The characterstic lengthscale L is de- 
fined by the minimum of the turbulent and Kolmogorov 
lengthscales 



L = C'l max 



^3/2 



e \ e 



1/4- 



with 



= 1.0 + LiniUi 



(15) 



(16) 



where rii is the unit wall-normal of the closest wall- 
element pointing outward of the fiow domain, while Cl 
and Cri arc model constants. The definition of in 
Eq. (16) signifies a slight departure from the original 
model by attributing anisotropic and wall-dependent be- 
havior to its value. In the case of a channel flow where 
the wall is aligned with x, the wall- normal n = (0, — 1, 0). 
This gives = 2.3 in the computation of P22 in Eq. (11) 
and = 1.0 for all other components. The modifica- 
tion improves the centerline behavior of the wall-normal 
Reynolds stress component (w^) and in turn the cross- 
stream mixing of the passive scalar. Another departure 
from the original model is the application of the ellip- 
tic term L^V^pij (as originally proposed by Durbin^) as 
opposed to LV^(Lpij). This simplification was adopted 
because no visible improvement has been found by em- 
ploying the second, numerically more expensive term. 

The right hand side of Eq. (11) may be any local model 
for pressure redistribution: here we follow Dreeben and 
Pope^^ and use the stochastic Lagrangian equivalent of 



a modified isotropization of production (IP) model pro- 
posed by Pope*^. Close to the wall, the elliptic term on 
the left hand side of Eq. (11) brings out the non-local, 
highly anisotropic behavior of the Reynolds stress tensor, 
whereas far from the wall the significance of the elliptic 
term vanishes and the local model on the right hand side 
is recovered. 

The above model needs to be augmented by an equa- 
tion for a quantity that provides length-, or time-scale 
information for the turbulence. With traditional moment 
closures the most common approach is to solve a model 
equation for the turbulent kinetic energy dissipation rate 
e itself as proposed by Hanjalic and Launder®. An al- 
ternative method is to solve an equation for the mean 



characteristic turbulent frequency 
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and to define 



(17) 



In PDF methods, however, a fully Lagrangian description 
has been preferred. A Lagrangian stochastic model has 
been developed for the instantaneous particle frequency 
uj by van Slooten et al."^ of which different forms exist, 
but the simplest formulation can be cast into 



dw = — C3 (uj) {uj — (uj)) dt — Suj (1^) L^dt 



1/2 



dW, 



(18) 



where S^^ is a source/sink term for the mean turbulent 
frequency 



Cull 



V 



(19) 



where V ~ — (uiUj) d (Ui) /dxj is the production of tur- 
bulent kinetic energy, dW is a scalar valued Wiener- 
process, while 6*3,(74, Ctji and C^2 are model constants. 
Since the no-slip condition would incorrectly force e to 
zero at the wall, Eq. (17) needs to be modified, thus the 
dissipation is defined as^"^ 



(c.) (fc + ,yC^ (io)) , 



(20) 



where Ct is also a model constant. A simplification of 
the original model for the turbulent frequency employed 
by Dreeben and Pope^'^ is the elimination of the ad-hoc 
source term involving an additional constant, since no 
obvious improvement has been found by including it. 

Similarly to the model equations for the Lagrangian 
velocity and frequency increments, the evolution equa- 
tions of the passive scalar are also given in Lagrangian 
form, which define the two micromixing models that are 
investigated in this study: 



d^j 
dtp 



1 



{tp-mdt 



(lEM), 



(21) 



tr 



(V'-(0|F))dt (lECM), (22) 



where tp represents the sample space variable of the 
species concentration (p, is the micromixing timescale. 
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while {(f>\V) = {(j)\U {x , t) = V) denotes the expected 
value of the mean concentration conditional on the ve- 
locity. Both of these models represent the physical pro- 
cess of dissipation and reflect the concept of relaxation 
towards a scalar mean with the characteristic timescale 
im- The difference is that in the lEM model, all particles 
that have similar position interact with each other, while 
in the lECM model only those particles interact that also 
have similar velocities, e.g. fluid elements that belong to 
the same eddy. 

It can be shown that in the case of homogeneous turbu- 
lent mixing with no mean scalar gradient the two models 
are equivalent*^ and the micromixing timescale tm is pro- 
portional to the Kolmogorov timescale t = k/e. In the 
inhomogeneous case of a concentrated source, however, 
there are various stages of the spreading of the plume 
requiring different characterizations of tm. In this case, 
the formal simplicity of the lEM and lECM models is a 
drawback, since a single scalar parameter tm has to ac- 
count for all the correct physics. The timescale should 
be inhomogeneous and should depend not only on the 
local turbulence characteristics but also on the source 
location, type, size, distribution and strength. Because 
of this complexity, a general flow- independent specifica- 
tion of t,„ has been elusive. In the following, we define a 
micromixing timescale for a passive scalar released from 
a concentrated source into an inhomogeneous flow sur- 
rounded by no-slip walls. 

In general, is assumed to be proportional to the 
timescale of the instantaneous plume'^*. Once the ini- 
tial conditions are forgotten, theoretical results^" show 
that the timescale of the instantaneous plume is linear 
in t in the inertial subrange and is proportional to the 
turbulence timescale in the far field, when the instanta- 
neous plume grows at the same rate as the mean plume. 
Based on these considerations the micromixing timescale 
is computed according to 



<m = min 



1/3 



Ct 



— , Ct 



where rg denotes the radius of the source, {U) ^ is the 
mean velocity at the centerline of the channel, while Cg 
and Ct are micromixing model constants. This definition 
rcfiects the three stages of the spreading of the plume. 
In the first stage, the time scale of the plume is pro- 
portional to that of the source^^: accordingly, the first 
term in the min operator represents the effect of the 
source. In the second stage im increases linearly as the 
scalar is dispersed downstream and the distance x from 
the source grows^°. In the final stage, the timescale is 
capped with the characteristic timescale of the turbu- 
lence, which provides an upper limit in the third term of 
Eq. (23). Following Durbin^^ this is defined as the max- 
imum of the turbulent and Kolmogorov timescales: far 
from the boundaries it becomes A:/e, whereas near a sur- 
face, where — > 0, the Kolmogorov timescale provides a 
lower bound as Ct{v / sY^"^ ■ 



The equations to model the joint PDF of velocity and 
turbulent frequency closely follow the work of Dreeben 
and Pope^^. Slight modifications consist of 

• the anisotropic definition of lengtliscale L in 
Eqs. (15) and (16), 

• the application of the elliptic term L'^S/'^pij instead 
of L\/^{Lpij) in Eq. (11), and 

• the elimination of an ad-hoc source term in 
Eq. (19). 

These modifications do not make the methodology less 
general or flow-dependent, i.e. the velocity model remains 
complete. On the other hand, the specification of in 
Eq. (23) is somewhat fiow-dependent, since we assume 
that the scalar plume will be dispersed in the x (down- 
stream) direction. 

In summary, the flow is represented by a large num- 
ber of Lagrangian particles representing a finite sample 
of all fluid particles in the domain. Each particle has po- 
sition Xi and carries its velocity Ui , turbulent frequency 
id and concentration tp. These particle properties are ad- 
vanced according to Eqs. (5), (8), (18) and either (21) 
or (22), respectively. The discretized particle equations 
are advanced in time by the first order accurate forward 
Euler method. This method was preferred to the more in- 
volved exponential scheme that was originally suggested 
by Dreeben and Pope^'^, since the code is sufficiently sta- 
ble with the simpler and computationally less expensive 
Euler method as well. 



III. NUMERICAL METHOD 

In turbulent channel flow after an initial development 
region, the statistics of fluid dynamics (e.g. turbulent 
velocity and frequency) are expected to become homo- 
geneous in the streamwise direction, while strongly in- 
homogeneous in the wall-normal direction. A passive 
scalar released into this flow from a concentrated source 
initially exhibits high inhomogeneity and only far down- 
stream becomes fully mixed. We will separately discuss 
the numerical issues related to the streamwise statisti- 
cally homogeneous one-dimensional fluid dynamics and 
the inhomogeneous two-dimensional scalar field. 



A. Modeling the fluid dynamics 

While the velocity and turbulent frequency statistics 
become one-dimensional, both the streamwise x and 
cross-stream y components of the particle position are 
retained so that particles can represent the streamwise 
inhomogeneity of the scalar. A one-dimensional grid, 
that is refined at the wall, is used to compute Eule- 
rian statistics of the velocity and turbulent frequency by 
ensemble averaging in elements. The elliptic-relaxation 
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equation (11) is also solved on this grid with a finite 
element method. A constant unit mean streamwise pres- 
sure gradient is imposed, drives the flow and builds up 
the numerical solution. The cross-stream mean pressure 
gradient is obtained by satisfying the cross-stream mean 
momentum equation 



l d(P) 
P dy 



dy 



(24) 



Wall-boundary conditions for the particles are the same 
as suggested by Dreeben and Pope^'^ but repeated here 
for clarity. Over any given time-interval a particle un- 
dergoing reflected Brownian motion in the vicinity of a 
wall may strike the wall infinitely many times^*^. There- 
fore wall-conditions have to be enforced on particles that 
either penetrate or possibly penetrate the wall during 
timcstep At, the probability of which can be calculated 
by53 



P = exp 



uAt J 



(25) 



where 3^q and y denote the distance of the particle from 
the wall at the previous and current timcstep, respec- 
tively. Thus particle wall-conditions are applied if 



or if 



y ^ and exp 



-y^y 



> V, 



(26) 



(27) 



where ij is a random variable with a standard uniform 
distribution. The imposed wall-conditions are 



Ui = (no slip) 



(28) 



and the particle frequency is sampled from a gamma dis- 
tribution with mean and variance 



1 dV2fc 
Ct dy 



and Ca 



(29) 



Only half of the channel is retained for the purpose of 
modeling the fluid dynamics. The particle conditions for 
the center line are symmetry conditions, i.e. if the par- 
ticle tries to leave the domain through the centerline, 
it is reflected back with opposite normal velocity. Con- 
sistently with the above particle conditions, boundary 
conditions are imposed on the Eulerian statistics as well. 
At the wall, the mean velocity and the Reynolds stress 
tensor are forced to zero. The mean frequency (w) is 
set according to Eq. (29). At the centerline, the shear 
Reynolds stress (uv) is set to zero. At the wall in the 
elliptic- relaxation equation (11), is set according to 
pij ~ —A.^eniUj. In the current case the wall is aligned 
with y ~ 0, thus only the wall-normal component is non- 
zero: p22 = — 4.5£. At the centerline, symmetry con- 
ditions are applied for p^, i.e. homogeneous Dirichlet- 
conditions are applied for the off-diagonal components 



and homogeneous Neumann-conditions for the diagonal 
components. 

An equal number of particles are uniformly generated 
into each element and initial particle velocities are sam- 
pled from a Gaussian distribution with zero mean and 
variance 2/3, i.e. the initial Reynolds stress tensor is 
isotropic with unit turbulent kinetic energy. Initial par- 
ticle frequencies are sampled from a gamma distribution 
with unit mean and variance 1/4. 



B. Modeling the passive scalar 

A passive, inert scalar is released from a concentrated 
source. The scalar statistics are inhomogeneous and in 
general not symmetric about the channel centerline, thus 
a second, two-dimensional grid is employed to calculate 
scalar statistics. The use of separate grids for the fluid 
dynamics and scalar fields enables the grid refinement 
to be concentrated on different parts of the domain, i.e. 
the scalar-grid can be refined around the source, while 
the fluid dynamics-grid is refined at the wall. The two- 
dimensional grid is unstructured and consists of triangles. 
The role of this mesh is twofold: it is used for comput- 
ing Eulerian scalar statistics and for tracking particles 
on the domain. Since the scalar is passive, only one- 
way coupling between the two grids is necessary. This is 
accomplished by using the local velocity statistics com- 
puted in the Id-elements for those 2d-elements which lie 
closest to them in the wall-normal coordinate direction. 
The Eulerian statistics in both grids are computed by 
a two-step procedure: first ensemble averaging is used 
to compute statistics in elements, then these element- 
based statistics are transferred to nodes by averaging the 
elements surrounding the nodes. When Eulerian statis- 
tics are needed in particle equations, the average of the 
nodal values are used for all the particles that reside in 
the given element. Spatial derivatives are computed with 
linear finite element shapefunctions for triangles^'*. The 
long two-dimensional rectangular domain is also subdi- 
vided into several equally sized bins. The velocity and 
turbulent frequency statistics are computed using the 
one-dimensional grid in which only particles in the first 
bin participate. The position of these particles are then 
copied to all downstream bins and mirrored to the upper 
half of the channel. If the particle hits the centerline, its 
concentration is exchanged with its mirrored pair on the 
upper half, allowing a possible non-symmetric behavior 
of the scalar. 

In order to numerically compute the expected value of 
the mean concentration conditional on the velocity field 
{4i\V) in Eq. (22), one needs to discretizc the velocity 
sample space V and compute different means for each 
sample space bin. A straightforward way to implement 
this is to equidistantly divide the three-dimensional ve- 
locity space into cubical bins and compute separate scalar 
means for each cube using those particles whose veloc- 
ity fall into the given cube. A better way of choosing 
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the conditioning intervals for each dimension is to define 
their endpoints so that the probabihties of the particle 
velocities falling into the bins are equal. For a Gaus- 
sian velocity PDF these endpoints can be obtained from 
statistical tables as suggested by Fox'*^ for homogeneous 
flows. If the approximate shape of the joint velocity PDF 
is not known, as in our case, a different sort of algorithm 
is required to homogenize the statistical error over the 
sample space. 

Note that the sample-spatial distribution of the con- 
ditioning intervals has to be neither equidistant nor the 
same in all three dimensions and can also vary from el- 
ement to clement. Therefore, wc determine the binning 
dynamically so that no bin remains without particles and 
the number of particles in each bin is approximately the 
same. A judicious sorting and grouping of the particles 
according to their velocity components can be used to 
achieve this. We compute conditioned means in each tri- 
angular element as follows. For a binning of 2 x 2 x 2 (i.e. 
the desired total number of bins is 8), first the particles 
residing in the triangle are sorted according to their lA 
velocity component. Then both the first and the second 
halves of the group are separately sorted according to 
their V velocity component. After further dividing both 
halves into halves again, the quarters are sorted accord- 
ing to the W component. Finally, halving the four quar- 
ters into eight subgroups, we compute scalar means for 
each of these eight subgroups. This procedure is general, 
it employs no assumptions about the shape of the veloc- 
ity PDF and homogenizes the statistical error over the 
sample space. It also ensures that there will be no empty 
bins as long as the number of conditioning bins and the 
number of particles/elements are reasonable. In princi- 
ple, the number of subdivisions (i.e. the sample-spatial 
refinement) can be arbitrary. We employ the binning 
structure of 4 x 4 x 4 throughout the current calculations. 

In order to reduce the statistical error during the com- 
putation of {4)\V), Fox^^ proposed an alternative method 
in which the three dimensional velocity space is projected 
onto a one dimensional sub-space where the sample- 
spatial discretization is carried out. This way, a relatively 
large number of particles can be used to obtain a local 
{(j)\V)^ which results in a smaller statistical error. The 
projection however is exact only in the case of a passive 
inert scalar in homogeneous turbulent shear flows with a 
uniform mean scalar gradient. In the current, more gen- 
eral inhomogeneous context, this projection could still 
be used as a modeling assumption''^. To explore the er- 
ror introduced by the projection in modeling inhomoge- 
neous flows, we implemented and tested this method for 
the channel flow. For both releases (to be discussed in 
the next sections) wc found that the projection has the 
largest influence on the scalar PDFs by homogenizing 
their high peaks, i.e. shaping them closer to a Gaussian. 
In the case of the centerline release, an artifact of double- 
peaks in the r.m.s. proflles of the cross-stream distribu- 
tion of the passive scalar can also be attributed solely to 
this projection method. In the following sections, only 



TABLE I: Constants for modeling the joint PDF of velocity 
and frequency. 
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C4 


Ct 


Cl 


Crj 




75 




1.85 


0.63 


5.0 


0.25 


6.0 


0.134 


72.0 


1.4 


0.1 


0.5 0.73 



the former, more general, three-dimensional method was 
used in conjunction with the lECM model with 4^^ bins. 



IV. RESULTS 

The model has been run for the case of fully devel- 
oped channel flow at Re-r = 1080 based on the friction 
velocity Ur and the channel half-width h with a pas- 
sive scalar released from a concentrated source at the 
centerline {us/h = 1.0) and in the viscous wall region 
{Us/h = 0.067). The results are divided into a discussion 
of the fluid dynamics statistics (IV A), a comparison of 
the two micromixing models (IV B) and a presentation of 
unconditional (IV C) and conditional (IV D) scalar statis- 
tics. 



A. Fluid dynamics 

The equations to model the velocity and turbulent fre- 
quency have been solved on a 100-cell one dimensional 
grid with 500 particles per cell. The applied model con- 
stants are displayed in Table I. The computed cross- 
stream profiles of mean streamwise velocity, the non-zero 
components of the Reynolds stress tensor and the rate 
of dissipation of turbulent kinetic energy are compared 
with the DNS data of Abe et al.^^ at Rcr = 1020 in 
Fig. 1. Previous PDF modeling studies employing the 
elliptic relaxation technique^^'^'^'^^ have been conducted 
up to Rct = 590. The high-level inhomogeneity and 
anisotropy in the viscous wall region are well represented 
by the technique at this higher Reynolds number as well. 
The purpose of including the parameter in Eq. (15) of 
the wall-normal component of pij is to correct the over- 
prediction of the wall-normal Reynolds stress component 
(w^) at the centerline. This facilitates the correct be- 
havior of the mean of the dispersed passive scalar in the 
center region of the channel (presented in Section IV B). 



B. Comparison of the lEM and lECM micromixing 
models 

An often raised criticism of the IBM model is that there 
is no physical basis for assuming the molecular mixing 
to be independent of the velocity field. This assump- 
tion gives rise to a spurious (and unphysical) source of 
scalar fiux^^. This behavior of the lEM model has also 
been demonstrated for line sources in homogeneous grid 
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FIG. 1: Cross-stream profiles of (a) tlie mean streamwise velocity, (b) the diagonal components of the Reynolds stress tensor, 
(c) the shear Reynolds stress and (d) the rate of dissipation of turbulent kinetic energy. Lines - PDF calculation, symbols - 
DNS data of Abe et al."^"^. All quantities are normalized by the friction velocity and the channel half- width. The DNS data is 
scaled from Rcr = 1020 to 1080. 



turbulence^^. The situation can be remedied by introduc- 
ing the velocity-conditioned scalar mean which 
leads to the lECM model. Often invoked as a desirable 
property of micromixing models is that the scalar PDF 
should tend to a Gaussian for homogeneous turbulent 
mixing^°'^^ (i.e. statistically homogeneous scalar field in 
homogeneous isotropic turbulence). While mathemati- 
cally a Gaussian docs not satisfy the boundcdness prop- 
erty of the advection-diffusion scalar transport equation, 
it is generally assumed that the limiting form of the PDF 
can be reasonably approximated by a clipped Gaussian. 
Also, Chatwin^'"'^® argued that in most practical cases, 
where the flow is inhomogeneous, the scalar PDF is better 
approximated by non-Gaussian functions, which should 
ultimately converge to a Dirac delta function about the 
mean, S(')p — (0)), where (0) approaches a positive value 
in bounded domains and zero in unbounded domains. 

In fully developed turbulent channel flow the center 
region of the channel may be considered approximately 
homogeneous"^^'^^. Thus for a center line source, up to a 
certain downstream distance where the plume still lies 
completely in the center region, the mean scalar field can 
be described by Taylor's theory of absolute dispersion^°. 



Likewise, numerical simulations are expected to repro- 
duce experimental measurements of grid turbulence. Ac- 
cording to the theory, the mean-square particle displace- 
ment (3^^) is related to the autocorrelation function of 
the Lagrangian velocity Rl = {v{t)v{t')) / (v^) as 



"'0 



(30) 



where it is assumed that in stationary turbulence Rj^ de- 
pends only on the time difference ^ = t — t' . Lagrangian 
statistics such as i?L(C) are difficult to determine experi- 
mentally. An analytical expression that is consistent with 
the theoretically predicted behavior of the Lagrangian 
spectrum in the inertial subrange is^^ 



i?i(0 = cxp(-M), 



(31) 



where denotes the Lagrangian integral timescale. 
Substituting Eq. (31) into Eq. (30) the following analyt- 
ical expression can be obtained for the root-mcan-square 
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FIG. 2: Cross-stream mean concentration profiles normalized by their respective peak values at difTerent downstream locations 
as computed by the (a) lECM and (b) lEM models for the centerline release. Lines - PDF calculation at solid line, x/h = 4.0, 
dashed line, x/h = 7.4 and dot-dashed line, x/h = 10.8, hollow symbols - analytical Gaussians using Eq. (34) at o, x/h = 4.0; 
A, x/h = 7.4 and □, x/h — 10.8, filled symbols - experimental data of Lavertu and Mydlarski'^^ at •, x/h = 4.0 and x/h = 7.4. 
Also shown, PDFs of scalar concentration fluctuations at {x/h — 7.4, y/h = 1.0) for the (c) lECM and (d) lEM models. Lines 
- computation, symbols - experimental data. 



particle displacement 



— 1 

Tl 



t 



(32) 



This expression can be used to approximate the spread of 
the plume that is released at the centerline of the channel. 
As Lagrangian timescale we take 



(33) 



where Co is usually taken as the Lagrangian veloc- 
ity structure function inertial subrange constant'^^^^^, 
which ensures consistency of the Langevin equation (8) 
with the Kolmogorov hypothesis in stationary isotropic 
turbulence^". In the current case the value of Co is de- 
fined by Eq. (10) and is no longer a constant, but depends 
on the velocity statistics. For the purpose of the cur- 
rent analytical approximation, however, a constant value 
(0.8) has been estimated as the spatial average of Cq com- 
puted by Eq. (10). For the cross-stream Reynolds stress 



(v^) and the dissipation rate e their respective centerline 
values are employed. In analogy with time t in homo- 
geneous turbulence, we define t = x/ (U)^, where x is 
the downstream distance from the source and (U)^ is the 
mean velocity at the centerline. Thus the cross-stream 
mean scalar profiles predicted by Eq. (30) are obtained 
from the Gaussian distribution 



(0(y)> 



Q 



1/2 



exp 



2al 



(34) 



where Q is the source strength and ys is the cross-stream 
location of the source. 

After the velocity field converged to a statistically sta- 
tionary state, a passive scalar is continuously released 
from a concentrated source. Two release cases have been 
investigated, where the scalar has been released at the 
centerline {ys/h = 1.0) and in the close vicinity of the 
wall {ys/h = 0.067). The viscous wall region experiences 
the most vigorous turbulent activity. The turbulent ki- 



10 



TABLE II: Model constants of the microniixing timescale t, 
defined by Eq. (23) for both the lEM and lECM models. 





source location 




Cs 


Ct 


centerline 


ys/h = 1.0 


y+ = 1080 


0.02 


0.7 


wall 


ys/h = 0.067 


y+ = 72 


1.5 


0.001 



netic energy, its production and its dissipation and the 
level of anisotropy all experience their peak values in this 
region, see also Fig. 1 (b). This suggests a significantly 
different level of turbulent mixing between the two re- 
lease cases. Accordingly, the constants that determine 
the behavior of the micromixing timescales have been se- 
lected differently. Both the lEM and lECM models have 
been investigated with the micromixing timescale defined 
by Eq. (23) using the model constants displayed in Table 
II. 

The different behavior of the two models is demon- 
strated in Fig. 2, which shows mean concentration pro- 
files for the centerline release computed by both the lEM 
and lECM models together with the analytical Gaus- 
sian solution (34) and the experimental data of Lavertu 
and Mydlarski'^^ for turbulent channel flow. Indeed, the 
downstream evolution of the cross-stream mean concen- 
tration profiles computed by the lECM model follows the 
Gaussians and is expected to deviate far downstream in 
the vicinity of the walls, where the effect of the walls is 
no longer negligible. It is also apparent in Fig. 2 (b) that 
the lEM model changes the mean concentration, as ex- 
pected. As discussed by Lavertu and Mydlarski'^-'^ , the 
measurements of the mean concentration experience the 
largest uncertainty due to inaccuracies in estimating the 
free-stream mean. Also, to improve the signal-to-noise 
ratio far downstream, a thicker wire had to be employed 
for measurements performed on the second half of the 
length considered, i.e. x/h > 11.0. These difficulties are 
probably the main source of the discrepancy between the 
experimental data and the agreeing analytical and nu- 
merical results for the case of the centerline release. Be- 
cause of these inconsistencies only results for the first 
half of the measured channel length (x/h < 11.0) are 
considered in the current study. 

The marginal PDF of scalar concentration can be ob- 
tained from the joint PDF f{V,ijj) by integrating over 
the velocity space 

fW = I f{V,ij)dV. (35) 

According to experimental data in grid turbulence'^® the 
skewness at the centerline is expected to be negative close 
to the source and to become positive only farther down- 
stream. At x/h = 7.4, y/h = 1.0 the temperature PDF 
measured by Lavertu and Mydlarski'^^ suggests positive 
skewness in accordance with Sawford's data"^®. In Fig. 2 
(c) and (d) the normalized PDFs of scalar concentration 
fluctuations at this location as computed by both models 
arc depicted. As opposed to the lEM model prediction, 



both the location of the peak and the overall shape of 
the PDF are captured correctly by the lECM model. 

The different behavior of the two micromixing models 
is apparent in all one point statistics considered, with the 
lECM model producing a closer agreement to experimen- 
tal data. The price to pay for the higher accuracy is an 
additional 30-40% in CPU time as compared to the lEM 
model. In the remaining section only the lECM model 
results are considered. 



C. Scalar statistics with the lECM model 

Cross-stream distributions of the first four moments 
of the scalar concentration at different downstream loca- 
tions are shown in Fig. 3 for both release scenarios. The 
results are compared to experimental data where avail- 
able. 

The mean and root-mean-square (r.m.s.) profiles are 
normalized by their respective peak values. The width of 
the mean concentration profiles is most affected by the 
wall-normal Reynolds stress component (v^) which is re- 
sponsible for cross-stream mixing. Due to the underpre- 
diction of this component by the velocity model through- 
out most of the inner layer {y~^ < 800) and the uncer- 
tainties in the experimental data mentioned in Section 
IV B, the mean concentration profiles in Fig. 3 should be 
considered at most qualitative. 

For the wall-release, the r.m.s. profiles display a clear 
drift of the peaks towards the centerline with increas- 
ing distance from the source Fig. 3 (f). This tendency 
has also been observed in turbulent boundary layers by 
Fackrell and Robins'^" and Raupach and Legg^^. Since 
the scalar is statistically symmetric, in the case of the 
centerline release, no tranverse drift of the r.m.s. pro- 
files is expected. Fig. 3 (b). Double peaking of the r.m.s. 
profiles has been observed in homogeneous turbulence by 
Warhaft^^ and Karnik and Tavoularis^®, noting that the 
profiles are initially double-peaked close to the source, 
then single-peaked for a short distance and then again 
double-peaked far downstream. Lavertu and Mydlarski'^^ 
found no double peaks in their measurements. Corre- 
sponding to the channel fiow experiments, the PDF sim- 
ulation exhibits no double-peaking in the r.m.s. profiles. 
Applying the projection method to compute {4>\V) as de- 
scribed in Section III B results in double peaking of the 
r.m.s. profiles, which is possibly due to a loss of statis- 
tical information due to its Gaussian assumption of the 
velocity PDF. 

Skewness profiles are depicted in Fig. 3 (c) and (g). 
For both release cases, near the centers of the plumes the 
skewness is close to zero, indicating that the PDFs of the 
scalar concentration downstream of the sources are ap- 
proximately symmetric. Towards the edges of the plumes 
however, the PDFs become very highly positively skewed, 
with a sudden drop to zero in the skewness outside of 
the plume. As observed by Lavertu and Mydlarski'^^, 
the downstream evolutions of the skewness profiles in- 
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FIG. 3: Cross-stream distributions of the first four moments of scalar concentration at different downstream locations for (a)- 
(d) the centerline release {ys/h = 1-0) and (e)-(h) the wall release {ys/h = 0.067). Lines - calculations, symbols - experimental 
data at solid line, •, x/h — 4.0; dashed line, x/h — 7.4 and dot-dashed line, ■, x/h — 10.8. The horizontal dashed lines 
for the skewness and kurtosis profiles indicate the Gaussian values of and 3, respectively. Note the logarithmic scale of the 
kurtosis profiles. 
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dicate the eventual mixing of the plume, with the high 
peaks decreasing. In the current simulations the high 
skewness-peaks at the edge of the plumes start increas- 
ing first to even higher levels (up to about x/h = 10.0) 
and only then start decreasing. In the case of the wall- 
release, the negative skewness in the viscous wall region 
(also apparent in the experimental data) becomes even 
more pronounced in the buffer layer and in the viscous 
sublayer, where experimental data is no longer available. 
The kurtosis values are close to the Gaussian value of 
3 at the cross-stream location of the sources, but show 
significant departures towards the edges of the plume. 



Figure 4 shows downstream evolutions of the peak of 
the mean and r.m.s. and the width of the mean con- 
centration profiles. In homogeneous isotropic turbulence 
and homogeneous turbulent shear flow the decay rate of 
the peak of the mean concentration profiles is reasonably 
well described by a power law of the form (0)p(,ak ^ 
In the present inhomogeneous fiow Lavertu and Myd- 
larski'^-^, based on the experiments, suggest decay expo- 
nents of n ~ —0.7 and —0.6 for the wall and centerline 
sources, respectively. These evolutions are compared to 
experimental data in Fig. 4 (a) and (b). Downstream 
evolutions of the width of the mean concentration pro- 
files Umcan are plotted in Fig. 4 (c) and (d) for the two 
releases. According to the experimental data, these do 
not exhibit power-law dependence, as is the case in ho- 
mogeneous flows. Since the simulations are carried out 
only on the first half of the measured channel length, 
the three downstream locations are not sufficient to un- 
ambiguously decide whether the simulation data exhibits 
power-law behavior for the peaks and widths of the mean 
profiles. 



The downstream decay of the peak values of the r.m.s. 
profiles can be well-approximated by a power-law of the 
form (^'^)p^(,^ak ^ ^ similarly to homogeneous shear flow 
and isotropic grid-generated turbulence. Fig. 4 (e) and 
(f). The experiments suggest n = — 1 for both releases. 



Probability density functions of scalar concentration 
fluctuations are depicted in Fig. 5 for both release cases. 
The cross-stream location of these PDFs are chosen to co- 
incide with that of their respective sources, i.e. y/h= 1.0 
for the centerline release and y/h = 0.067 for the wall- 
release. Two downstream locations are plotted, at the 
first and at the third location from the sources mea- 
sured by Lavertu and Mydlarski'^^, at x/h = 4.0 and 
x/h = 10.8, respectively. While the PDFs for the cen- 
terline release are in reasonable agreement with the ex- 
periments, some discrepancies are apparent in the wall- 
release case. A possible reason behind this disparity is the 
ad- hoc specification of the mixing timescale in Eq. (23), 
which is mostly based on theoretical considerations and 
experimental observations in homogeneous turbulence. 



D. Conditional statistics 

The current model solves for the full joint PDF of 
the turbulent velocity, frequency and scalar concentra- 
tion. Therefore we can also examine those quantities that 
require closure assumptions in composition-only PDF 
methods. This methods are often used in combustion 
engineering to model complex chemical reactions in a 
given turbulent flow or in dispersion modeling in the at- 
mospheric boundary layer. In these cases the simplest 
approach is to assume the shape of the velocity PDF and 
numerically solve a set of coupled model equations that 
govern the evolution of the joint PDF of the individual 
species concentrations in composition space. The conser- 
vation equation for a single reactive scalar is 



(36) 



where S{(j)) is the chemical source term. In high Reynolds 
number, constant property flow the PDF of a reactive 
scalar f^{ip;x,t) is governed by^°''^^ 



di 
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(37) 
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An attractive feature of these formulations is that the 
usually highly nonlinear chemical source term is in closed 
form. Closure assumptions, however, are necessary for 
the velocity fluctuations conditional on the scalar con- 
centration {ui\ip) and the conditional scalar dissipa- 
tion (2rV0 • V(/)|'(/') or the conditional scalar diffusion 
(rV^^lV')- Since for the current case S{(f>) ~ 0, the 
marginal scalar PDF /(V') defined in Eq. (35) is equal to 
f(f,, thus in the remaining text we use them interchange- 
ably. 

For the convective term Dopazo^'* applied the linear 
approximation 



(39) 



to compute the centerline evolution of the temperature 
PDF in a turbulent axisymmetric heated jet. This lin- 
ear approximation is exact for joint Gaussian velocity 
and scalar fiuctuations. While many experiments^^"^^ 
confirm the linearity of the conditional mean veloc- 
ity around the local mean conserved scalar, Kuznetsov 
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FIG. 4: Downstream evolutions of (a), (b) the peak mean scalar concentration, (c), (d) the width of the mean concentration 
and (e), (f) the peak of the r.m.s. profiles for the centerline and wall releases, respectively. Solid lines - numerical results, 
symbols - experimental data. 



and Sabel'nikov''^ observe that most of the experimental 
data show departure from this linear relationship when 
\i) — (0)1 is large. Experimental data from Sreenivasan 
and Antonia^^ and Bilger et al.^^ also show that in inho- 
mogcncous flows the joint PDF of velocity and scalar is 
not Gaussian, which makes the above linear approxima- 
tion dubious in a general case. Nevertheless, this linear 
model is sometime applied to inhomogeneous scalar fields 
because of its simplicity. 

Another commonly employed approximation is to in- 



voke the gradient diffusion hypothesis 

-/^(^,|^) = rT|^, (40) 

where Tt{x^ t) is the turbulent diffusivity. In the current 
case, we specify the turbulent viscosity vt based on the 
traditional k — e closure and relate it to F^ with the 
turbulent Prandtl number ctt as 

r. = ^ = (41) 
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FIG. 5: Probability density functions of scalar concentration fluctuations at selected downstream locations for the (a) centerline 
and (b) wall-releases at the cross-stream location of their respective sources (i.e. y/h — 1.0 and y/h — 0.067, respectively). 
Lines - calculation, symbols - experimental data at solid line, •, x/h = 4.0 and dot-dashed line, ■, x/h = 10.8. 




FIG. 6: Cross-stream velocity fluctuation conditioned on the scalar concentration for the wall-release {ys/h = 0.067). Thick 
lines, lECM model; thin lines, gradient diffusion approximation of Eq. (40); straight sloping lines, linear approximation of 
Eq. (39). (a) downstream evolution at the height of the source y/h — 0.067: solid lines, x/h = 4.0; dot-dashed lines, x/h = 10.8 
and (b) cross-stream evolution at x/h = 7.4: solid lines, y/h = 0.067; dot-dashed lines, y/h = 0.67. 



where C^i = 0.09 is the usual constant in the fc — e model 
and we choose ctt = 0.8. 

In Fig. 6 (a) the downstream evolution of the cross- 
stream velocity fluctuation conditioned on the scalar is 
depicted for the wall-release case. Both locations are at 
the height of the source, i.e. y/h = 0.067. The concen- 
tration axis for both locations is scaled between their 
respective local minimum and maximum concentration 
values, Tpmin and T/imax- Note that the model curves show 
higher negative velocity for low-concentration particles 
as the distance from the source increases. This is ex- 
pected, since particles deep inside the plume can have 
very low concentrations only if they did not come from 
the source but traveled very fast from above, so that they 
did not have much time to exchange concentration with 
the source material. As the plume spreads, only parti- 
cles with stronger negative velocity can maintain their 
low concentration values. Likewise, as the center of the 



plume moves towards the centerline of the channel, high- 
concentration particles also need to have stronger nega- 
tive velocities to escape from exchange during their jour- 
ney from the plume-center to our sensors, which is ap- 
parent on the right side of the figure. Obviously, the 
linear approximation (39) cannot be expected to cap- 
ture the non-linearity of the model curves, but except 
for extremely low and high concentrations it performs 
reasonably well. On the other hand, the gradient diffu- 
sion approximation is capable of capturing most features 
of the lECM model behavior: it successfully reproduces 
the non-linearity, with some discrepancy at low and high 
concentrations. It is also apparent that the numerical 
computation of the derivatives of the PDFs in the gra- 
dient diffusion model (40) is most sensitive to sampling 
errors at the concentration extremes due to lower number 
of particles falling into the concentration bins there. 

The cross-stream evolution of the conditioned velocity 
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FIG. 7: lECM model predictions for the mean scalar dissipation conditioned on the concentration for (a) the centerline release 
(jjs/h = 1.0) and (b) the wall-release {ys/h = 0.067) at different downstream locations: solid line, x/h = 4.0; dashed line, 
x/h = 7.4 and dot-dashed line, x/h = 10.8. The cross-stream locations are the same as the respective source positions. Note the 
different scales for the dissipation curves between the different releases. Also shown are the scalar PDFs at the same locations 
for both releases in (c) and (d), respectively. 



fluctuation is shown in Fig. 6 (b). Both sensors are now 
at the downstream location x/h — 7.4 with increasing 
distance from the wall at y/h = 0.067 and 0.67. As the 
sensor moves towards the channel centerline, the detected 
low-concentration particles need weaker negative velocity 
to maintain those low concentrations. The sensor loca- 
tions relative to the plume centerline can be identified by 
examining the cross-stream velocity of the high concen- 
tration particles. The sensors at y/h = 0.067 and 0.67 
are below and above the plume centerline, respectively, 
since high-concentration particles at these locations pos- 
sess negative and positive cross-stream velocities. As 
is expected, the linear approximation reasonably repre- 
sents the model behavior for mid-concentrations, while 
its performance degrades at locations with higher non- 
Gaussianity, i.e. towards the edge of the plume. The 
performance of the gradient diffusion model is reason- 
able, except at the concentrations extremes. 

For the lECM micromixing model, the mean dissipa- 
tion conditioned on the scalar concentration can be com- 
puted from^^ 



where 



0(7A)= / {(b\V)f{V\^)dV. 



(43) 



The function (/)('0) in Eq. (43) can be obtained by 
taking the average of {4>\V) over those particles that 
reside in the bin centered on ip. The integral in 
Eq. (42), however, is more problematic. As Sawford^^ 
notes, numerical integration errors that accumulate at 
extreme concentrations may be amplified when divided 
by the scalar PDF approaching zero at those locations. 
Since the integral over all concentrations vanishes, i.e. 
(2r(V0)2|-i/'max) /^(V'max) = 0, for mid-conccntrations it 
can be evaluated either from the left (V'min — ^ V') or from 
the right ('0max — >■ Thus the integration errors at the 
concentration extremes can be significantly decreased by 
dividing the domain into two parts, integrating the left 
side from the left and the right side from the right and 
merging the two results in the division-point. Due to 
statistical errors, however, the integral over all concen- 
trations may not vanish. In that case, the nonzero value 
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FIG. 8: Mean scalar difTusion conditioned on the concentration as predicted by the lECM and lEM models for (a) the centerline 
release {ys/h = 1.0) and (b) the wall-release {t/s/h = 0.067) at different downstream locations. The cross-stream locations are 
the same as the respective source positions. Solid line, x/h = 4.0; dashed line, x/h = 7.4 and dot-dashed line, x/h = 10.8. The 
straight lines are the linear predictions of the lEM model of Eq. (46). 



can be distributed over the sample space by correcting 
the integrand with the appropriate fraction of this error 
in each bin. 

The conditional mean dissipation for three different 
downstream locations is depicted in Fig. 7 for both re- 
lease cases. As for the conditional velocity, the abscissas 
here are also scaled between the local V'min and V'max- 
The dissipation is normalized by the mixing timescale 
and the square of the mean scalar peak ((/>)pga^i;. at 
the corresponding downstream locations. Note that in 
the case of the wall-release, the dissipation curves are an 
order of magnitude lower than in the centerline release 
case. This is mainly a result of the choice of the different 
micromixing model constants, especially Ct- 

In the case of the wall release, the curves exhibit bi- 
modal shapes at all three downstream locations. This 
tendency has also been observed by Kailasnath ct al.^'^ 
in the wake of a cylinder and by Sawford in a double- 
scalar mixing layer^^ and, to a lesser extent, also in ho- 
mogeneous turbulence''^. Sardi et al.''''^ suggest that in 
assumed-PDF methods of turbulent combustion a qual- 
itative representation of the conditional dissipation can 
be obtained in terms of the inverse PDF. To examine this 
relationship, the corresponding scalar PDFs arc also plot- 
ted in Fig. 7 with the same scaling on the concentration 
axis as the dissipation curves. It is apparent that these 
results support this reciprocal connection except at the 
extremes: high values of the PDF correspond to low dis- 
sipation (and vice versa). This can be observed for both 
releases, but it is most visible in the wall-release case, 
where the mid-concentration minimum between the two 
maxima of the bi-modal dissipation curves correspond to 
the peaks in the PDFs. 

The lECM model (22) implies a model for the mean 
diffusion conditioned on the scalar concentration as 

(rv20|v,V') = --^(^-(0|y)). (45) 



The downstream evolution of the conditional diffusion is 
depicted in Fig. 8 for both releases. The concentration 
axes are scaled as before and the curves are normalized 
by the scalar variance (0'^), the concentration at the 
source, 0o and the mean unconditioned dissipation (x) = 
(2r(V(?!))^), which is computed by integrating Eq. (42) 
over the whole concentration space. Also shown are the 
predictions according to the lEM model, which is given 
by the linear relationship"^^ 




Far downstream as the scalar gets better mixed, the pre- 
dictions of the IBM and lECM models get closer. This 
behavior has been observed for other statistics, as well as 
for other flows such as the double-scalar mixing layer"^^. 
Kailasnath et al.'''^ report experimental data on similar 
shapes for the conditional diffusion in the turbulent wake 
of a cylinder. 

V. CONCLUSIONS 

Several different stochastic models have been com- 
bined to develop a complete PDF-IECM model to com- 
pute the joint PDF of turbulent velocity, frequency and 
scalar concentration in fully developed turbulent chan- 
nel flow. The flow is represented by a large number of 
Lagrangian particles and the governing stochastic differ- 
ential equations have been integrated in time in a Monte- 
Carlo fashion. The high anisotropy and inhomogeneity 
at the low-Reynolds-number wall-region have been cap- 
tured through the elliptic relaxation technique, explic- 
itly modeling the vicinity of the wall down to the viscous 
sublayer by imposing only the no-slip condition. Durbin^ 
suggested the simple LRR-IP closure of Launder et al.^, 
originally developed in the Eulerian framework, as a local 
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model used in the elliptic relaxation equation (11). Since 
then, several more sophisticated local Reynolds stress 
models have been investigated in conjunction with the 
elliptic relaxation technique''^. In the PDF framework, 
the Lagrangian modified IP model of Pope'*^ is based on 
the LRR-IP closure. We introduced an additional model 
constant in the definition of the characteristic length- 
scale L (15) whose curvature determines the behavior of 
the relaxation and, ultimately, the overall performance of 
the model in representing the Reynolds stress anisotropy. 
This resulted in a correction of the original model over- 
prediction of the wall- normal component (w^) far from the 
wall, which crucially influences the cross-stream mixing 
of the transported scalar. However, increasing the con- 
stant adversely affects the level of anisotropy that can 
be represented by the technique. A more accurate treat- 
ment of the Reynolds stresses and scalar mixing should 
be achieved by a more elaborate second moment closure, 
such as the nonlinear C-L model of Craft and Launder*"^ 
or the Lagrangian version of the SSG model of Speziale 
et al.^ suggested by Pope^^. 

An unstructured triangular grid is used to compute 
Eulerian scalar statistics and to track particles through- 
out the domain. The main purpose of employing unstruc- 
tured grids has been to prepare the methodology for more 
complex flow geometries. A similar particle-in-cell ap- 
proach has been developed by Muradoglu et al.^'''''^ and 
by Zhang and Haworth^^ for the computation of turbu- 
lent reactive flows. These approaches combine the ad- 
vantages of traditional Eulerian CFD codes with PDF 
methods in a hybrid manner. Our aim here is to develop 
a method that is not a hybrid one, so the consistency 
between the computed fields can be naturally ensured. 
The emphasis is placed on generality, employing numeri- 
cal techniques that assume as little as possible about the 
shape of the numerically computed joint PDF. 

We compared the performance of the lEM and the 
lECM micromixing models in an inhomogeneous flow 
with strong viscous effects by modeling both the tur- 



bulent velocity field and the scalar mixing. The more 
sophisticated lECM model provides a closer agreement 
with experimental data in channel flow for the additional 
computational expense of 30-40% compared to the lEM 
model. 

Several conditional statistics that often require closure 
assumptions in PDF models where the velocity field is 
assumed were extracted and compared to some of their 
closures. In particular, our conclusions suggest that the 
scalar-conditioned velocity is well approximated by a lin- 
ear assumption for mid-concentrations at locations where 
the velocity PDF is moderately skewed. The gradient dif- 
fusion approximation, however, captures most features 
including the nonlinearity and achieves a closer agree- 
ment with the lECM model in slightly more skewed re- 
gions of the flow as well. At local concentration extremes 
and in extremely skewed regions the gradient diffusion 
approximation markedly departs from the lECM model. 
The mean scalar dissipation conditioned on the scalar 
concentration may be well-approximated by the inverse 
relationship suggested by Sardi et al.''^ in inhomogeneous 
flows with significant viscous effects as well, except at the 
concentration extremes. In computing the conditional 
scalar diffusion, both the lEM and the lECM models 
produce similar slopes due to the same scalar dissipation 
rate attained. 



Acknowledgments 

We wish to thank Thomas Dreeben for the helpful and 
insightful discussions on the velocity model, Rodney Fox 
for his interesting comments on the lECM model, Hi- 
royuki Abe and Hiroshi Kawamura for providing the DNS 
data on the turbulent velocity field and Laurent Myd- 
larski, who kindly provided the experimental data on the 
scalar statistics. 



* Electronic address: jbakosi@gmu.edu 

^ D. C. Haworth and S. B. Pope, "A generalized Langevin 
model for turbulent flows", Phys. Fluids 29, 387 (1986), 
URL http : //link . aip . org/link/?PFL/29/387/l . 

^ P. A. Durbin, "A Reynolds stress model for near- wall tur- 
bulence", J. Fluid Mech. 249, 465 (1993). 

^ P. R. van Slooten, Jayesh, and S. B. Pope, "Ad- 
vances in PDF modeling for inhomogeneous turbu- 
lent flows", Phys. Fluids 10, 246 (1998), URL 
http : //link . aip . org/link/?PHF/10/246/ 1 . 

* W. P. Jones and B. E. Launder, "The prediction of lami- 
narization with a two-equation model of turbulence". Int. 
J. Heat Mass Tran. 15, 301 (1972). 

^ D. P. Bacon, N. N. Ahmad, Z. Boybeyi, T. J. Dunn, M. S. 
Hall, P. C. S. Lee, R. A. Sarma, M. D. Turner, K. T. W. 
HI., S. H. Young, et al., "A dynamically adapting weather 
and dispersion model: The Operational Multiscale Envi- 



ronment Model with Grid Adaptivity (OMEGA)", Mon. 
Weather Rev. 128, 2044 (2000). 
^ J. C. Rotta, "Statistiche theorie nichthomogener turbu- 
lenz", Z. Phys. 129, 547 (1951). 

B. E. Launder, G. J. Reece, and W. Rodi, "Progress in 

the development of a Reynolds-stress turbulent closure", 

J. Fluid Mech. 68, 537 (1975). 
^ K. Hanjalic and B. E. Launder, "A Reynolds stress model 

of turbulence and its application to thin shear flows", J. 

Fluid Mech. 52, 609 (1972). 
^ C. G. Speziale, S. Sarkar, and T. B. Gatski, "Modelling 

the pressure-strain correlation of turbulence: an invariant 

dynamical systems approach", J. Fluid Mech. 227, 245 

(1991). 

" S. B. Pope, Turbulent flows (Cambridge University Press, 

Cambridge, 2000). 
^ R. O. Fox, Computational models for turbulent reacting 



18 



flows (Cambridge University Press, 2003). 

S. B. Pope, "Computations of turbulent combustion: 

progress and challenges", Proc. Combust. Inst. 23, 591 

(1990). 

J. C. R. Hunt and J. M. R. Graham, "Free-stream tur- 
bulence near plane boundaries", J. Fluid Mech. 84, 209 
(1978). 

E. R. van Driest, "On the turbulent flow near a waU", J. 
Aeronaut. Sci. 23, 1007 (1956). 

Y. G. Lai and R. M. C. So, "On near wall turbulent flow 
modelling", J. Fluid Mech. 221, 641 (1990). 
^® T. J. Craft and B. E. Launder, "Reynolds stress closure 
designed for complex geometries". Int. J. Heat Fluid Fl. 
17, 245 (1996). 

W. Rodi and N. N. Mansour, "Low Reynolds number k — 
e modelling with the aid of direct numerical simulation 
data", J. Fluid Mech. 250, 509 (1993). 
B. E. Launder and D. B. Spalding, "The numerical com- 
putation of turbulent flows", Comput. Meth. Appl. M. 3, 
269 (1974). 

^® A. K. Singhal and D. B. Spalding, "Predictions of two- 
dimensional boundary layers with the aid of the k — e model 
of turbulence" , Comput. Meth. Appl. M. 25, 365 (1981). 
W. Rodi, Turbulence Models and Their Application in Hy- 
draulics - A State-of-the-art Review (Intl. Assoc. for Hy- 
draulics Research, Delft, The Netherlands, 1980). 
D. B. Spalding, GENMIX - A General Computer Program 
for Two-Dimensional Parabolic Phenomena (Intl. Assoc. 
for Hydraulics Research, 1977). 

T. D. Dreeben and S. B. Pope, "Probability den- 
sity function and Reynolds-stress modeling of near-wall 
turbulent flows", Phys. Fluids 9, 154 (1997), URL 
http : //link . aip . org/link/?PHF/9/154/l . 
T. D. Dreeben and S. B. Pope, "Probability density func- 
tion/Monte Carlo simulation of near- wall turbulent flows", 
J. Fluid Mech. 357, 141 (1998). 

B. I. Shraiman and E. D. Siggia, "Scalar turbulence". Na- 
ture 405, 639 (2000). 
'^^ Z. Warhaft, "Passive scalars in turbulent flows", Annu. 
Rev. Fluid Mech. 32, 203 (2000). 

Z. Warhaft, "The interference of thermal flelds from line 
sources in grid turbulence", J. Fluid Mech. 144, 363 
(1984). 

H. Stapountzis, B. L. Sawford, J. C. R. Hunt, and R. E. 
Britter, "Structure of the temperature fleld downwind of 
a line source in grid turbulence", J. Fluid Mech. 165, 401 
(1986). 

B. L. Sawford, "A simple representation of a developing 
contaminant concentration fleld", J. Fluid Mech. 289, 141 
(1995). 

U. Karnik and S. Tavoularis, "Measurements of heat diffu- 
sion from a continuous line source in a uniformly sheared 
turbulent flow", J. Fluid Mech. 202, 233 (1989). 
J. E. Fackrell and A. G. Robins, "Concentration fluctua- 
tions and fluxes in plumes from point sources in a turbulent 
boundary layer", J. Fluid Mech. 117, 1 (1982). 
R. A. Lavertu and L. Mydlarski, "Scalar mixing from a 
concentrated source in turbulent channel flow", J. Fluid 
Mech. 528, 135 (2005). 

R. D. Moser, J. Kim, and N. N. Mansour, "Di- 
rect numerical simulation of turbulent channel flow up 
to Rcr = 590", Phys. Fluids 11, 943 (1999), URL 
http : //link . aip . org/link/?PHF/ll/943/l . 
H. Abe, H. Kawamura, and Y. Matsuto, "Surface heat-flux 



fluctuations in a turbulent channel flow up to Re-r = 1020 
with Pr = 0.025 and 0.71", Int. J. Heat Fluid Fl. 25, 404 
(2004). 

A. J. Vrieling and F. T. M. Nieuwstadt, "Turbulent disper- 
sion from nearby point sources - interference of the con- 
centration statistics", Atmos. Environ. 37, 4493 (2003). 

J. Villermaux and J. C. Devillon, Representation de 
la coalescence et de la redispersion des domaines de 
segregation dans un fluide par un modele d'interaction 
phenomenologique, in Proceedings of the 2nd Intl. Symp. 
on Chemical Reaction Engineering (Elsevier, New York, 
1972), pp. 1-13. 

C. Dopazo and E. E. O'Brien, "An approach to the au- 
toignition of a turbulent mixture" , Acta Astronaut. 1, 1239 
(1974). 

C. Dopazo, Recent developments in pdf methods, in Turbu- 
lent reactive flows, edited by P. A. Libby (Academic, New 
York, 1994), pp. 375-474. 

B. L. Sawford, "Micro-mixing modelling of scalar fluctua- 
tions for plumes in homogeneous turbulence". Flow Tur- 
bul. Combust. 72, 133 (2004). 

B. L. Sawford, "Lagrangian modeling of scalar 
statistics in a double scalar mixing layer", Phys. 
Fluids 18, 085108 (pages 11) (2006), URL 
http : //link . aip . org/link/?PHF/18/085108/l . 
A. K. Luhar and B. L. Sawford, "Micromixing modelling of 
mean and fluctuating scalar flelds in the convective bound- 
ary layer", Atmos. Environ. 39, 6673 (2005). 
M. Cassiani, P. Franzese, and U. Giostra, "A PDF mi- 
cromixing model of dispersion for atmospheric flow. Part 
I: development of the model, application to homogeneous 
turbulence and to neutral boundary layer", Atmos. Envi- 
ron. 39, 1457 (2005). 

■'^ M. Cassiani, P. Franzese, and U. Giostra, "A PDF mi- 
cromixing model of dispersion for atmospheric flow. Part 
II: application to convective boundary layer" , Atmos. En- 
viron. 39, 1471 (2005). 

■'^ M. Cassiani, A. Radicchi, and J. D. Albertson, "Mod- 
elling of concentration fluctuations in canopy turbulence" , 
Boundary-Layer Meteorol. 122, 655 (2007). 
S. B. Pope, "PDF methods for turbulent reactive flows". 
Prog. Energ. Combust. 11, 119 (1985). 

'^^ C. W. Gardiner, Handbook of stochastic methods for 
physics, chemistry and the natural sciences (Springer- 
Verlag, Berlin Heidelberg New York, 2004), 3rd ed. 

'^^ A. Einstein, Investigations on the theory of the brownian 
movement (Methuen, London, 1926), translation by A. D. 
Cowper. 

S. B. Pope, "On the relationship between stochas- 
tic Lagrangian models of turbulence and second- 
moment closures", Phys. Fluids 6, 973 (1994), URL 
http : //link . aip . org/link/?PHF/6/973/l . 

D. C. Wilcox, Turbulence modeling for CFD (DCW Indus- 
tries, La Canada, CA, 1993). 

R. O. Fox, "On velocity-conditioned scalar mixing in ho- 
mogeneous turbulence", Phys. Fluids 8, 2678 (1996), URL 
http : //link . aip . org/link/?PHF/8/2678/l . 
P. Franzese and M. Cassiani, "A statistical theory of tur- 
bulent relative dispersion", J. Fluid Mech. 571, 391 (2007). 
G. K. Batchelor, "Diffusion in a fleld of homogeneous tur- 
bulence II. The relative motion of particles", P. Camb. 
Philos. Soc. 48, 345 (1952). 

P. A. Durbin, "Near-wall turbulence modeling without 
damping functions" , Theor. Comp. Fluid Dyn. 3, 1 (1991). 



19 



I. Karatzas and S. E. Shreve, Brownian motion and 

stochastic calculus (Springer, 1991), 2nd ed. 

R. Lohner, Applied CFD techniques: an introduction based 

on finite element methods (John Wiley & Sons. U.K., 

2001). 

M. Waclawczyk, J. Pozorski, and J.-P. Minier, "Probabil- 
ity density function computation of turbulent flows with a 
new near- wall model", Phys. Fluids 16, 1410 (2004), URL 
http : //link . aip . org/link/?PHF/16/1410/l . 

^® S. B. Pope, "The vanishing effect of molecular diffusivity 
on turbulent dispersion: implications for turbulent mixing 
and the scalar flux", J. Fluid Mech. 359, 299 (1998). 

^'^ P. C. Chatwin, "Some remarks on modelling the PDF of 
the concentration of a dispersing scalar in turbulence", 
Eur. J. Appl. Math. 13, 95 (2002). 

P. C. Chatwin, "Singular PDFs of a dispersing scalar in 
turbulence". Flow Turbul. Combust. 72, 273 (2004). 

■''^ G. Brethouwer and F. T. M. Nieuwstadt, "DNS of mixing 
and reaction of two species in a turbulent channel flow: a 
validation of the conditional moment closure" , Flow Tur- 
bul. Combust. 66, 209 (2001). 

®° G. I. Taylor, "Diffusion by continuous movements", P. 
Loud. Math. Soc. 20, 196 (1921). 

P. S. Arya, Air pollution meteorology and dispersion (Ox- 
ford University Press, 1999). 

®^ A. S. Monin and A. M. Yaglom, Statistical fluid mechanics, 
vol. 2 (The MIT Press, 1975). 

®^ M. R. Raupach and B. J. Legg, "Turbulent dispersion 
from an elevated line source: measurements of wind- 
concentration moments and budgets", J. Fluid Mech. 136, 
111 (1983). 

®* C. Dopazo, "Probability density function approach 
for a turbulent axisymmetric heated jet. Center- 
line evolution", Phys. Fluids 18, 397 (1975), URL 
http : //link . aip . org/link/?PFL/18/397/ 1 . 
V. A. Bezuglov, Measurement of velocity of concentra- 
tion transport, Tech. Rep. Pt. 1, Dolgoprudnyi, Moscow 
Physical-Technical Institute (1974). 

®^ Y. V. Golovanov, Experimental investigation of statisti- 
cal characteristics of turbulent fluctuations of contaminant 
concentration in a submerged jet, Ph.D. thesis, Dolgoprud- 
nyi, Moscow Physical- Technical Institute (1977). 

®^ Y. A. Shcherbina, Statistical characteristics of turbulent 



transport. Tech. Rep., Dolgoprudnyi, Moscow Physical- 
Technical Institute (1982). 

K. S. Venkataramani and R. Chevray, "Statistical features 
of heat transfer in grid generated turbulence: constant gra- 
dient case", J. Fluid Mech. 86, 513 (1978). 
K. R. Sreenivasan and R. A. Antonia, "Joint probability 
densities and quadrant contributions in a heated turbulent 
round jet", AIAA J. 16, 867 (1978). 
™ V. R. Kuznetsov and V. A. Sabel'nikov, Turbulence and 
combustion (Hemisphere, New York, 1990), English ed. 
R. W. Bilger, L. R. S. L. V., and Krishnamoorthy, "Re- 
action in a scalar mixing layer", J. Fluid Mech. 233, 211 
(1991). 

B. L. Sawford, "Conditional scalar mixing statistics in ho- 
mogeneous isotropic turbulence". New J. Phys. 6, 1 (2004). 
P. Kailasnath, K. R. Sreenivasan, and J. R. Saylor, "Con- 
ditional scalar dissipation rates in turbulent wakes, jets, 
and boundary layers", Phys. Fluids 5, 3207 (1993), URL 
http : //link . aip . org/link/?PFA/5/3207/l . 
■^"^ K. Sardi, A. M. K. P. Taylor, and J. H. Whitelaw, "Condi- 
tional scalar dissipation statistics in a turbulent counter- 
flow", J. Fluid Mech. 361, 1 (1998). 

V. Whizman, D. Laurence, M. Kanniche, P. A. Durbin, 
and A. Demuren, "Modeling near-wall effects in second- 
moment closures by elliptic relaxation". Int. J. Heat Fluid 
Fl. 17, 255 (1996). 
^® T. J. Craft and B. E. Launder, Computation of impinging 
flows using second moment closures, in Proc. 8th sympo- 
sium on Turbulent Shear Flows, T. U. Munich, Germany 
(1991), pp. 5-8. 

M. Muradoglu, P. Jenny, S. B. Pope, and D. A. Caughey, 
"A consistent hybrid finite- volume/particle method for the 
PDF equations of turbulent reactive fiows", J. Comput. 
Phys. 154, 342 (1999). 

M. Muradoglu, S. B. Pope, and D. A. Caughey, "The hy- 
brid method for the PDF equations of turbulent reactive 
fiows: consistency conditions and correction algorithms", 
J. Comput. Phys. 172, 841 (2001). 

Y. Z. Zhang and D. C. Haworth, "A general mass con- 
sistency algorithm for hybrid particle/finite- volume PDF 
methods", J. Comput. Phys. 194, 156 (2004). 



